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With  attractive  advantages  of  high  efficiency,  energy  saving  and  environmental  friendliness,  the  ground 
source  heat  pump  (GSHP)  system  has  been  used  widely  in  China  in  recent  years.  This  paper  summarizes 
the  analytical  solution,  numerical  solution  and  experimental  investigation  of  the  heat  transfer  of  the 
ground  heat  exchanger  (GHE),  analyzes  the  simulation  model  and  long-term  operation  performance  of 
the  GSHP  system,  and  introduces  the  lastest  hybrid  ground  source  heat  pump  (HGSHP)  system.  In 
addition,  this  paper  discusses  and  summarizes  the  shortages  and  imperfects  of  the  current  research  on 
the  simulation  of  the  GSHP  system  and  gives  some  recommendations  for  future  work. 
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1.  Introduction 

As  the  global  energy  crises  and  environmental  problems  become 
more  and  more  serious,  the  ground  source  heat  pump  (GSHP) 
technology  has  showed  a  trend  of  booming  growth  in  world  due 
to  its  attractive  advantages  of  high  efficiency,  energy  saving  and 
environmental  friendliness  and  the  GSHP  system  has  become  a 
hotspot  in  clean  energy  research.  In  recent  years,  great  numbers  of 
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the  GSHP  systems  have  been  used  in  commercial  and  residential 
buildings  throughout  China.  According  to  the  Ministry  of  Science 
and  Technology  of  China,  the  total  area  adopting  the  GSHP  systems 
will  be  up  to  350  million  square  meters  by  the  end  of  2015  [1].  Now 
in  China,  the  GSHP  technology  is  in  the  stage  of  rapid  development. 


2.  Heat  transfer  model  and  operation  performance  of  GHE 

Ground  heat  exchanger  (GHE)  is  an  important  part  of  the  GSHP 
system  [2].  The  heat  transfer  of  the  GHE  is  a  very  complicated 
dynamic  process.  On  one  hand,  the  heat  transfer  of  the  GHE 
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Nomenclature 

a 

thermal  diffusivity  (m2/s) 

p 

I  integration  variable  (-) 

q 

heat  flux  (W/m2) 

p 

density  (kg/m3) 

T 

temperature  (K) 

Cp 

specific  heat  at  constant  pressure  (J/((kg  K)) 

Subscripts 

r 

distance  from  the  borehole  axis  (m) 

t 

time  (s) 

0 

initial 

m 

mass  flow  rate  of  the  circulating  fluid  in  the 

/ 

fluid 

U-tube  (kg/s) 

5 

solid 

R 

heat  transfer  resistant  ((m  1<)/W) 

p 

ground  heat  exchanger 

L 

borehole  length  (m) 

b 

borehole  wall 

u 

velocity  of  circulating  fluid  in  the  U-tube  (m/s) 

in 

inlet  of  the  U-tube 

Fo 

Fourier  constant  (-) 

out 

outlet  of  the  U-tube 

a 

air 

Greek  letters 

inlet  branch  of  the  U-tube 

2 

outlet  branch  of  the  U-tube 

1 

thermal  conductivity  (W/(m  K)) 

depends  on  the  buried  method,  soil  characteristics  underground 
hydrological  parameters,  backfill  materials  and  meteorological 
data.  On  the  other  hand,  the  heat  transfer  of  the  GHE  interplays 
with  the  operation  performance  of  heat  pump  unit  and  the  load  of 
buildings  [3],  GHEs  are  divided  into  two  categories:  horizontal 
GHEs  and  vertical  GHEs.  Since  the  horizontal  GHE  covers  larger 
area  and  the  heat  transfer  process  is  considerably  influenced  by 
the  ground  surface  temperature  and  ambient  air  temperature, 
this  paper  focuses  on  the  vertical  GHE. 


method  to  obtain  the  temperature  response  of  the  finite  line  heat 
source  in  the  semi-infinite  medium  and  the  temperature  response 
is  expressed  as: 


Ts-T 


f 


[erfc 

y/p2+(Z~hf 

2y/a.st 

\  ^ 

s2+(z-h)2 

erfc 


Jpj+jz+h)2 


2s/ast 


\Jps2  +  (z+h)2 


\dh 


(2) 


where  z  is  the  axis  of  the  vertical  direction,  h  is  the  integration 
variable  of  the  finite  line  heat  source. 


2.1.  Analytical  solution 


2.1.2.  Cylindrical  source  theory 


2.1.1.  Line  source  theory 

In  1948,  Ingersoll  et  al.  [4]  improved  Kelvin’s  line  source 
theory  to  solve  the  heat  transfer  problem  of  the  GHE.  In  Ingersoll’s 
theory,  the  heat  transfer  of  the  GHE  is  simplified  as  the  heat 
transfer  of  the  heat  source  which  has  the  same  axis  of  the 
borehole.  The  model  has  several  assumptions  listed  as  follows: 
the  ground  is  an  infinite  medium  with  the  initial  uniform 
temperature  and  its  thermophysical  properties  are  homogeneous 
and  stable:  the  heat  transfer  in  the  direction  of  the  borehole  axis, 
including  the  heat  flux  across  the  ground  surface  and  down  the 
bottom  of  the  borehole,  is  neglected;  the  geometric  dimensions  of 
the  borehole  are  neglected  and  assumed  as  a  line  source  of  the 
borehole  axis. 

According  to  the  Ingersoll’s  line-source  theory,  the  temperature 
response  in  the  ground  due  to  a  constant  heat  flux  is  given  by: 


Ts-T0  = 


(1) 


where  I(X)  is  the  exponential  integral  function,  X  =  r/(2v/o^t);  Ts  is 
the  ground  temperature,  K;  T0  is  the  initial  ground  temperature,  K\ 
q  is  the  constant  heat  flux  of  the  line-source,  W/m;  r  is  the  distance 
from  the  borehole  axis,  m;  As  is  the  thermal  conductivity  of  the 
ground,  W/(m  K),  as  is  the  thermal  diffusivity  of  the  ground,  m2/s; 
t  is  the  time  since  the  beginning  of  operation,  s;  [1  is  integration 
variable. 

The  investigation  indicated  that  the  third  assumption  is  not 
reasonable  when  the  time  is  very  short  (t<  5rj;/as,  several  hours 
in  general).  For  the  GHE  with  constant  thermal  capacity,  the 
analytical  solution  shows  major  errors,  especially  for  short  time. 

The  traditional  infinite  line  source  model  ignored  the  influence 
of  the  finite  length  of  the  borehole  and  the  boundary  condition  of 
the  ground  surface.  Fang  et  al.  [5]  utilized  the  Green  function 


2. 1.2.1.  Cylindrical  source  theory  with  a  constant  heat  flux.  Cylindrical 
source  theory  extends  a  line  source  to  a  cylindrical  source  with  a 
constant  radius.  Since  the  analytical  solution  of  the  cylindrical 
source  model  has  distinct  physical  meanings,  and  the  cylindrical 
source  theory  is  more  accurate  than  the  line  source  theory,  the 
cylindrical  source  theoiy  is  more  popular  in  application  and  plays  as 
the  foundation  of  the  great  majority  of  numerical  simulation 
models. 

In  1954,  Ingersoll  et  al.  [6]  developed  the  analytical  solution  of 
the  cylindrical  source  with  a  constant  heat  flux  and  the  solution  is 
described  as: 


ATg  —  Tw—Tg  — 


q  G(Fo,P) 

L  J.s 


(3) 


where  q  is  the  heat  flux  of  the  GHE,  W;  L  is  the  borehole  length,  m; 
Fo  is  Fourier  constant;  P  is  the  ratio  between  the  distance  from 
the  borehole  axis  and  the  borehole  radius;  2S  is  the  mean  thermal 
conductivity  of  the  ground,  m2/s.  As  defined  by  Carslaw  et  al.  [7], 
the  expression  G  (Fo,  p)  is  only  a  function  of  time  and  distance 
from  the  borehole  center. 


2. 1.2.2.  Cylindrical  source  theory  with  varied  heat  fluxes.  In  actual 
operation,  the  heat  extracted  or  released  by  the  GHE  varies  as  the 
operation  conditions.  Eq.  (1)  should  be  improved  since  it  is 
derived  on  the  base  of  the  constant  heat  flux.  Bernier  [8] 
utilized  the  step  heat  flux  to  resolve  the  problem  of  the  varied 
heat  flux.  Any  load  which  changes  with  the  time  is  regarded  as 
the  superposition  of  the  thermal  effect  on  the  borehole  caused  by 
several  piecewise  linear  step  heat  fluxes.  Thus,  the  solution  of 
temperature  distribution  is  treated  as  the  temperature  response 
of  an  infinite  medium  caused  by  a  series  of  step  load  with 
different  heat  fluxes. 
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For  the  time  tn,  the  temperature  difference  between  far 
boundary  soil  and  borehole  wall  is  given  as: 


A7g,  =  IT £  £  <7 i{G[Fo(t,-ti-\)\  — G[Fo(tj— fj)] } 


i  =  ij  =  i 


(4) 


At  the  time  point  of  tn,  the  temperature  depends  on  both  the 
heat  flux  of  the  current  time  point,  and  the  heat  flux  of  previous 
time  points.  The  investigation  indicates  that  the  heat  flux  of  the 
current  time  point  plays  a  major  role  in  the  current  temperature 
while  the  effect  of  the  heat  flux  of  previous  time  points  is  much 
smaller.  And  the  further  the  previous  time  point  away  from  the 
calculation  time  point  is,  the  smaller  the  impact  is.  Thus,  calculat¬ 
ing  the  temperature  and  heat  flux  of  previous  time  points  one 
after  the  other  is  the  basis  of  obtaining  the  temperature  and  heat 
flux  of  the  current  time  point.  Since  it  is  difficult  to  calculate  the 
performance  characteristics  of  the  GHE  during  a  whole  year  or  a 
longer  period,  Yavuzturk  et  al.  [9]  introduced  the  term  of  the  load 
aggregation  to  deal  with  the  effect  of  historical  heat  flux  on  the 
current  temperature  distribution.  Given  that  the  effect  of  the 
current  heat  flux  is  larger  than  that  of  the  previous  heat  flux, 
during  simulation,  from  the  time  point  of  tm,  the  heat  flux  is  the 
superposition  of  the  heat  flux  of  the  subsequent  time  points. 
While  for  the  time  point  before  tm,  the  heat  flux  is  considered  as 
the  mean  heat  flux  of  the  previous  time  points,  qmean.m<  and  the 
superposition  is  not  necessary.  The  mean  heat  flux,  qmean,m.  is 
named  as  the  load  aggregation  and  the  period  from  tm  to  tn  is 
called  as  the  dominant  period  of  the  heat  influence.  Therefore,  the 
temperature  of  the  time  point  of  tn  can  be  expressed  as: 

ATg.i  =  Tg,i-Tw,i  =  jF  |qmean,i-a{G[fo(t,-t0)-G[Fo(ti-ti_/1)]}  +  ■  ■  ■ 

+  J2q‘~A+i  x  { G[F o(Gi— tj_i )— G[F o(t,A— tj)] } 
j'  =  i 

where  A  is  the  number  of  terms  of  the  superposition.  Eq.  (5) 
shows  that  the  number  of  terms  of  the  superposition  in  the 
temperature  calculation  equation  at  the  time  point  t„  reduces 
A-l  by  introducing  the  definition  of  the  load  aggregation.  The 
temperature  mainly  depends  on  the  dominant  time  of  the  heat 
influence. 

The  cylinder  source  theories  with  a  constant  heat  flux  and 
variable  heat  fluxes  both  assume  that  for  the  same  time  and  the 
same  operating  mode,  the  heat  flux  of  buried  tubes  with  different 
lengths  is  the  same  and  the  heat  flux  of  the  left  and  right  side 
tubes  with  the  same  depth  is  the  same,  which  is  inconsistent  with 
the  practical  condition.  In  the  computational  domain  coupling  the 
heat  transfer  of  the  buried  tubes  and  soil,  the  temperature 
gradient  of  the  buried  tube  wall  is  the  largest.  Therefore,  the 
assumption  of  the  cylinder  source  theory  can  produce  a  certain 
amount  of  errors.  Yavuzturk  et  al.  [10]  simplified  the  boundary 
condition  of  the  buried  tube  wall  as  a  constant  heat  flux  and 
defined  that  the  heat  flux  in  the  left  side  tube  is  40%  of  the  total 
heat  flux  and  the  right  side  is  60%.  Although  this  method  has 
some  degree  of  randomness,  it  is  reasonable  in  a  certain  extent 
because  it  indicates  the  heat  flux  inhomogeneity  between  the  left 
tube  and  right  tube. 

Among  domestic  research  institutions  focusing  on  the  analy¬ 
tical  solution  of  the  heat  transfer  of  the  GHE,  the  main  one  is 
Shandong  Jianzhu  University.  Researchers  of  the  institution 
divided  the  heat  transfer  area  into  two  parts:  inside  borehole 
and  outside  the  borehole.  The  heat  transfer  process  of  the  soil 
outside  the  borehole  is  treated  as  the  transient  heat  transfer  of 
the  line  source.  Literatures  [11-14],  respectively  analyzed  both 
the  infinite  line  source  model  and  finite  line  source  model.  The 


heat  transfer  of  the  inner  of  the  borehole,  including  backfill 
materials,  borehole  wall  and  circulating  fluid  inside  borehole,  is 
regarded  as  a  steady  heat  transfer  process  [15-19].  On  the  basis  of 
the  energy  equilibrium  theory,  literature  [15]  adopted  the  linear 
superposition  principle  to  establish  the  two-dimensional  heat 
conduction  model  during  the  flowing  process  of  fluid,  based  on 
which  the  dimensionless  temperature  of  the  fluid  along  the 
dimensionless  depth  of  the  borehole  in  vertical  U-tube  GHEs 
was  obtained.  Literature  [16]  considered  the  temperature  differ¬ 
ence  of  the  fluid  along  the  depth  of  the  borehole,  established  the 
quasi  three-dimensional  heat  transfer  model  of  the  inner  of  the 
borehole  and  derived  the  analytical  solution  of  the  vertical  GHE. 
Fig.  1  shows  the  heat  transfer  network  of  the  borehole  in  the  quasi 
three-dimensional  heat  transfer  model.  Fang  et  al.  [17]  also 
analyzed  the  heat  transfer  of  the  borehole  inner  of  the  double 
U-tube  GHE.  He  utilized  the  superposition  principle  to  establish 
the  quasi  three-dimensional  heat  transfer  model  of  the  borehole 
inner  and  derive  the  analytical  solution  of  the  vertical  double 
U-tube  GHE.  Given  that  in  the  heat  transfer  model  of  the  borehole 
inner,  the  borehole  wall  temperature  is  one  of  the  most  important 
parameters,  he  suggested  that  the  integration  mean  temperature 
along  the  whole  borehole  wall  in  the  depth  direction  can  be 
treated  as  the  representative  temperature  of  the  borehole  wall. 
Literature  [19]  adopted  the  virtual  heat  source  method  and  the 
linear  superposition  principle  to  obtain  the  steady  analytical 
solution  of  the  inclined  GHE. 

The  quasi  three-dimensional  heat  transfer  model  of  GHE  is 
expressed  as: 

I  mcnf^afi  =  <Ji+<Ji2  =  ik(Tb-Tfi(z))+j^(Tf2(z)-Tfi(z)) 

\  ~mcPfAIssr1=(l  2-<Ji2  =  wJTh-Tf2(z))--^-2(Tf2  (Z)  -Tfl  (Z)) 

Boundary  condition: 

Tf\(0)  =  Tfin,  Tfi  ( L)  =  Tf2(L )  (7) 

where  cpf  is  the  specific  heat  of  the  circulating  fluid  at  constant 
pressure,  J/( kg  C),  m  is  the  mass  flow  rate  of  the  circulating  fluid 
in  the  U-tube,  kg/s;  Tf  i  and  Tf2  are  the  inlet  temperature  and 
outlet  temperature  of  the  circulating  fluid  at  the  depth  of  z, 
respectively, °C;  Tb  is  the  borehole  wall  temperature,  /(;  Ri,  R2  and 
R12  are  the  heat  transfer  resistance  between  the  fluid  and  the 
borehole  wall  and  the  equivalent  heat  transfer  resistance  between 
adjacent  branch  pipes,  respectively,  R]  =R2,  (m  K)/W. 

By  means  of  Laplace  transform,  the  dimensionless  solution  of 
the  fluid  temperature  distribution  along  the  length  direction  of 
two  branch  pipes  are  derived  and  described  as: 

el(2)=/1(Z)+/2(Z)^G±Mi 

02(Z)  =  (  ’ 


borehole  wall(Tb  ) 


Fig.  1.  Heat  transfer  network  of  borehole 
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where 


_  z  ^  Tf](z)-Tb  ^  _  Tf2(z)-Tb 


Tfin~T, 


R\  = 

y 


rilpfR  1 


b 

mcpfR2 


Tfin~T[ 


b 

mcpfRu 


1 

(R*) 


'*'2  +  R*R*2 

1 


f1(Z)  =  ch(yZ)~ 


1 

y  \R*+R*2 


/2(Z)=ish(yZ) 

f3(Z)  =  ch(yZ)+l 


1 


y\R*+R*2J 


)sh(yZ) 


■)sh(yZ) 


Yang  et  al.  [20]  utilized  the  cylindrical  source  theory  with 
varied  heat  fluxes  to  solve  the  temperature  distribution  of  the 
borehole  wall.  Wang  et  al.  [21]  adopted  the  cylindrical  source 
theory  with  varied  heat  fluxes  to  establish  the  coupled  model 
combining  the  heat  pump  unit  and  the  GHE  and  simulated  the 
performance  characteristics  during  a  short  period  of  time. 

Above  studies  regarded  the  heat  transfer  of  the  GHE  as  a  pure 
heat  conduction  problem  and  neglected  the  effect  of  the  advec- 
tion  of  the  groundwater.  The  wet  soil  containing  water  has  two 
kinds  of  influence:  one  is  the  static  effect  of  the  humidity;  the 
other  is  the  effect  of  infiltration  or  flow  caused  by  the  humidity 
gradient.  The  humidity  in  the  ground  soil  converts  the  simple 
heat  conduction  problem  to  a  complicated  problem  concerning 
heat  convection  and  the  diffusion  of  heat  and  humidity  [22].  Diao 
et  al.  [23,24]  assumed  that  the  ground  was  the  uniform 
porous  medium  and  the  heat  transfer  of  the  ground  included 
the  heat  conduction  of  the  medium  and  the  heat  convection  of  the 
fluid.  The  advection  of  fluid  in  the  porous  medium  is  regarded  as  a 
moving  heat  source  and  the  analytical  solution  of  the 
two-dimensional  temperature  distribution  caused  by  the  heat 
source  in  the  infinite  medium  with  advection  was  resolved. 


2.2.  Numerical  solution 


Compared  with  the  analytical  solution,  although  the  numerical 
solution  cannot  show  the  distinct  physical  meaning  and  the  effect 
of  each  factor,  it  can  resolve  more  complicated  heat  transfer 
problems  of  the  GHE.  And  that’s  why  a  great  number  of  research¬ 
ers  simulated  numerically  the  heat  transfer  process  of  the  GHE. 
Literatures  [25-29]  adopted  an  infinite  cylinder  to  be  the  equiva¬ 
lent  of  the  U-tube  GHE  and  established  the  coupled  thermal 
transfer  differential  equations  of  the  GHE.  Fig.  2  shows  the  heat 
transfer  network  of  the  borehole.  The  corresponding  mathema¬ 
tical  expressions  are  as  follows: 


dt 


+  u 


57> 

8z 


gf 

PfCrf 


x<ra 


(9) 


dTp 

dt 


r  dr 


dr  dz 


dz 


ra<x<rb 


(10) 


_  dTs  1  d  , ,  dTs  8  , ,  dTs  n 
P'C„s—  =  )+  —  (lj— -)  0  <z<L,rb<r  <ra 

rs  p  dt  r  dr  dr  dz  dz  0 


Initial  condition: 


Ts(r,z,  0)  =  Tf(z,  0)  =  Tp(r,z,  0)  =  T0 
Fluid  inlet  temperature: 

Tf(z  =  0,t)  =  Tfin 


(ID 


(12) 

(13) 


Fig.  2.  Heat  transfer  network  of  borehole 


Boundary  condition: 


•0A| 
dr  I  r  =  r  . 


sni 

az  lz  =  i 


ha(Ta~Ts(z  —  0))  —  2S  |z  _  o 


(14) 


where  u  is  the  velocity  of  circulating  fluid,  m/s;  ra  and  rb  are  the 
equivalent  internal  and  external  radius  of  the  ground  heat 
exchanger,  respectively,  m;  is  the  length  of  far  boundary,  m. 
Subscript  /,  p,  and  s,  respectively  stand  for  fluid,  ground  heat 
exchanger  and  soil. 

Tang  et  al.  [30,31]  utilized  the  finite  volume  method  to 
simulate  numerically  the  temperature  and  flow  distribution  of 
the  single  U-tube  GHE.  In  the  axial  direction,  the  one-dimensional 
heat  convection  model  is  adopted  and  the  circular  tube  is 
replaced  by  the  fan-shaped  tube  in  order  to  obtain  convenient 
gridding.  In  the  depth  direction,  the  two-dimensional  heat  trans¬ 
fer  model  is  used  to  calculate  the  coupled  heat  transfer  perfor¬ 
mance  between  the  liquid  inside  the  borehole  and  the  soil.  The 
heat  transfer  models  of  these  two  directions  are  combined  into 
the  quasi  three-dimensional  heat  transfer  model.  During  the 
simulation,  the  heat  transfer  in  soil  is  assumed  to  be  conduction 
only  in  the  horizontal  direction;  the  effect  of  the  bend  tube  at  the 
bottom  of  U-tube  pipes  is  ignored,  and  the  fluid  at  the  bottom  of 
two  branch  pipes  is  assumed  to  have  the  same  temperature  and 
opposite  direction  of  velocity.  Zhao  et  al.  [32]  also  adopted  the 
one-dimensional  heat  transfer  model  in  the  vertical  direction  and 
the  two-dimensional  heat  transfer  model  in  the  horizontal  direc¬ 
tion.  Unlike  Tang’s  simulation  method,  he  developed  the  heat 
transfer  equations  of  two  leg  pipes  in  the  bipolar  coordinate 
system  respectively  and  calculated  the  actual  temperature  by 
overlapping  the  temperature  distributions  of  the  two  leg  pipes. 
On  the  previous  basis  and  in  consideration  of  the  interact  effect 
between  heat  and  humidity  transfer,  literature  [33]  set  up  the 
coupled  heat  and  humidity  transfer  model  of  the  soil  around  the 
U-tube  vertical  GHE  in  bipolar  cylindrical  coordinates.  Ren  et  al. 
[34]  extended  the  transient  two-dimensional  heat  transfer  model 
of  the  single  vertical  U-tube  GHE  proposed  by  Rottmayer  [35] 
and  analyzed  the  effect  of  the  heat  transfer  in  the  vertical 
direction  and  the  ground  surface  temperature  on  the  heat  transfer 
of  the  GHE.  Wang  [36]  established  a  mathematic  model  for  a  three- 
dimensional  heat  transfer  temperature  field  of  the  vertical  GHE 
in  the  ground  source  heat  pump  system.  He  presented  the  layer 
heat  transfer  theory,  in  which  the  vertical  GHE  and  the  surround¬ 
ing  ground  are  divided  into  three  heat  transfer  layers,  including 
the  saturated  heat  transfer  layer,  unsaturated  heat  transfer  layer 
and  none  heat  transfer  layer.  Yang  et  al.  [37]  developed  a  heat 
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exchange  model  for  the  vertical  U-shaped  buried  tubes  on  the  base 
of  the  element  energy  balance  method.  The  model  had  taken 
account  of  the  liquid  temperature  change  along  the  flow  path,  and 
reflects  thermal  interference  between  two  tube  feet  through  the 
introduction  of  a  thermal  interference  angle  and  an  equivalent  heat 
transfer  spacing.  Yang  et  al.  [38-41]  used  both  the  steady  and 
transient  heat  transfer  method  to  analyze  the  heat  transfer  char¬ 
acteristics  inside  and  outside  the  borehole  by  the  cylindrical  source 
with  variable  heat  fluxes.  Then  he  developed  a  heat  transfer  model 
of  the  soil  around  GHEs  by  taking  soil  freezing  into  consideration 
and  solved  the  freezing  phase  change  problem  by  the  apparent  heat 
capacity  method.  He  et  al.  [42]  treated  the  steady  heat  transfer  of 
the  soil  temperature  field  as  a  discontinuous  media  problem  with 
mixed  boundary  values.  Based  on  the  nonoverlapping  and  parallel 
domain  decomposition  method  of  boundary  element,  the  computa¬ 
tional  domain  is  divided  into  two  subdomains.  In  the  interface  of  the 
two  subdomains,  the  temperature  and  heat  flux  are  continuous. 

In  addition,  some  researchers  combined  the  analytical  solution 
and  numerical  solution  to  obtain  a  simple  and  quick  computa¬ 
tional  process.  Lu  et  al.  [43]  divided  the  computational  domain 
into  four  parts:  liquid  inside  the  borehole,  borehole  wall,  backfill 
materials  and  soil  of  the  numerical  computational  domain  and 
nonnumerical  computational  domain.  The  temperature  distribu¬ 
tion  of  the  liquid  inside  borehole,  borehole  wall,  backfill  materials 
and  soil  of  numerical  computational  domain  are  calculated  by  the 
numerical  iterative  method,  while  the  soil  temperature  distribu¬ 
tion  of  the  nonnumerical  computational  domain  is  calculated  by 
the  finite-length  linear  source  method.  The  boundary  tempera¬ 
ture  of  the  numerical  computational  domain  is  determined  by  the 
finite-length  linear  source  method,  while  the  heat  flux  per  unit 
length  is  obtained  by  calculating  the  convective  heat  transfer 
between  the  liquid  inside  the  borehole  and  the  borehole  wall. 
Thus  the  calculations  of  both  inside  and  outside  the  numerical 
computational  domain  are  coupled. 

Yuan  et  al.  [44]  divided  the  computational  domain  into  two 
parts  based  on  the  boundary  of  the  borehole  wall.  The  computa¬ 
tional  method  is  the  steady  analytic  solution  for  the  inner  of  the 
borehole  and  the  finite  volume  method  for  the  soil.  The  simula¬ 
tion  models  of  the  two  parts  are  coupled  via  the  borehole  wall 
temperature  and  the  heat  flux  through  the  borehole  wall.  By 
means  of  Laplace  transform,  the  outlet  temperature  of  the 
circulating  fluid  and  the  heat  per  unit  length  of  the  soil  can  be 
derived.  The  computational  process  is  divided  into  three  steps. 
First,  calculate  the  outlet  temperature  of  the  circulating  fluid  and 
heat  per  unit  length  based  on  the  initial  borehole  wall  tempera¬ 
ture.  Second,  calculate  the  temperature  distribution  of  the  soil 
under  the  variable  heat  flux  boundary  conditions,  then  calculate 
the  borehole  wall  temperature,  outlet  temperature  of  the  circu¬ 
lating  fluid  and  heat  per  unit  length  at  the  current  time.  Third,  use 
the  heat  flux  of  the  soil  at  the  current  time  as  the  heat  flux  of  the 


soil  for  the  next  time  and  then  to  repeat  the  previous  computa¬ 
tional  process. 

Fan  et  al.  [28,45]  analyzed  the  numerical  solution  of  the  GHE 
with  groundwater  advection.  Under  the  assumption  that  the 
velocity  of  the  groundwater  was  even  and  the  groundwater  only 
flowed  along  the  horizontal  direction,  the  simulation  results 
indicated  that  groundwater  flow  had  significant  influence  on  the 
original  temperature  contour,  and  the  larger  the  groundwater 
velocity  was,  the  more  significant  the  effect  was. 

In  practical  projects,  the  GHE  is  always  made  of  multiple  pipes. 
Different  from  the  heat  transfer  of  the  single  pipe,  the  heat 
transfer  of  each  pipe  in  multiple  pipes  influences  each  other. 
The  heat  transfer  of  multiple  pipes  is  very  complicate.  According 
to  reported  literatures,  the  heat  transfer  of  multiple  pipes  mainly 
focuses  on  the  temperature  response  of  the  soil.  Most  researchers 
adopted  the  constant  heat  flux  cylindrical  source  model  [46,47] 
and  discussed  the  effect  of  different  kinds  of  arrangements  on  the 
temperature  distribution  of  the  ground  soil  [48,49].  In  addition, 
Ding  et  al.  [50]  defined  the  correction  factor  of  multiple  pipes, 
which  is  the  ratio  of  the  heat  transfer  resistance  between  multiple 
pipes  and  the  ground  to  the  heat  transfer  resistance  between  the 
single  pipe  and  the  ground.  Then  he  calculated  the  correction 
factors  of  multiple  pipes  under  different  operation  modes. 

Due  to  the  complexity  of  the  geometry  of  the  vertical  U-tube 
GHE  and  the  coupled  heat  transfer  between  the  fluid  inside  pipes 
and  the  soil,  it  is  impossible  and  also  unnecessary  to  establish  an 
accurate  model  to  simulate  the  practical  condition  to  obtain  the 
accurate  solution.  Hence,  all  the  simulations  need  a  certain 
simplification.  Table  1  summarizes  the  assumed  conditions  of 
parts  of  the  referred  literatures.  However,  there  are  few  reports 
about  the  errors  of  the  temperature  distribution  inside  and 
outside  pipes  when  the  simplifications  are  used  to  solve  the 
solution. 

2.3.  Experimental  investigation 

Experimental  investigation  is  one  of  the  important  methods 
for  investigating  the  GSHP  system.  Wang  et  al.  [51]  carried  out 
comparative  experiment  for  the  coaxial  and  U-tube  GHE.  Experi¬ 
mental  results  showed  that  the  heat  transfer  performance  of  the 
coaxial  GHE  was  better  than  that  of  the  U-tube  type.  Zhao  et  al. 
[52]  indicated  that  the  better  the  insulation  of  the  inner  pipe  was, 
the  larger  the  heat  extraction  of  the  GHE  was.  Experimental 
results  also  proved  that  it  was  feasible  to  use  beads  glass  ball  to 
imitate  the  convective  soil  layer  around  the  GHE  Gao  et  al.  [53] 
analyzed  a  new  GHE  with  the  screw  core  tube  bundle.  The 
experimental  results  showed  that  the  heat  transfer  capacity  of 
the  new  GHE  was  as  high  as  two  times  of  that  of  the  U-tube  GHE. 
Li  et  al.  [54]  conducted  experiments  on  the  heat  extraction 
and  rejection  of  the  single  U-tube  and  double  U-tube  vertical 


Table  1 

Assumed  conditions  of  literatures. 


Assumed  conditions 

No.  of  literatures 

9 

10 

11 

12 

13 

16 

17 

18 

19 

20 

21 

24 

Ignoring  the  heat  conduction  transfer  of  soil  in  the  vertical  direction 

X 

X 

X 

X 

X 

X 

X 

V 

X 

V 

X 

X 

Same  flow  rate  of  circulating  fluid  in  any  cross-section 

V 

X 

X 

X 

X 

X 

V 

X 

X 

V 

X 

X 

Uniform  initial  temperature  of  ground 

X 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

Ignoring  the  thermal  effect  of  groundwater  advection 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

Ignoring  the  influence  of  ambient  air  temperature  fluctuations 

V 

V 

V 

V 

V 

V 

V 

X 

V 

V 

V 

V 

Cylindrical  heat  source 

X 

X 

V 

V 

V 

V 

V 

X 

V 

V 

X 

V 

Constant  thermal  properties  of  soil 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

Two-dimensional  heat  transfer 

V 

X 

X 

X 

V 

V 

V 

X 

X 

V 

X 

V 

Steady  heat  transfer 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

V 

X 
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GHE.  The  results  showed  that  single  U-tube  vertical  GHE  was 
better  than  double  U-tube  vertical  GHE  in  heat  extraction  and 
rejection.  Zhao  et  al.  [55]  tested  the  U-tube  pile  GHE  and  U-tube 
well  GHE  in  summer  and  found  that  the  heat  extraction  rate  of 
the  U-tube  pile  GHE  was  as  large  as  116%  of  that  of  the  U-tube 
well  GHE.  At  the  mean  time,  the  effective  thermal  radius  of  the 
U-tube  pile  GHE  was  also  larger  than  that  of  the  U-tube  well  GHE. 
Wang  et  al.  [56]  carried  out  comparative  investigation  on  the 
single  HDPE  U-shaped  tube,  single  HDPE  &  stainless  steel  U- 
shaped  tube  and  double  HDPE  U-shaped  tube.  Xue  et  al.  [57] 
indicated  that  the  heat  flux  of  the  double  U-shaped  tube  is  more 
than  that  of  the  single  U-shaped  tube  under  both  the  heat 
extraction  and  rejection  conditions.  Zhao  et  al.  [58]  established 
an  in-situ  field  thermal  response  test  on  the  basis  of  the 
cylindrical  heat  source  model. 


3.  Simulation  model  and  operation  performance  of  GSHP 

3.J.  Coupled  model  and  simulation  of  GHE  and  heat  pump  unit 

During  the  operation  of  the  GSHP  system,  the  practical  opera¬ 
tion  performance  of  the  GHE  directly  affects  the  performance  of 
the  heat  pump  unit.  Thus,  in  the  GSHP  system,  the  dynamic  load 
of  the  building,  the  performance  of  the  heat  pump  unit  and  the 
GHEs  influence  each  other. 

In  China,  in  terms  of  the  simulation  of  the  GSHP  system,  the 
studies  focus  on  the  simulation  of  the  heat  pump  unit  or  a  single 
component,  and  the  coupling  of  the  building  load,  ground  heat 
pump  unit  and  GHE.  Since  the  operating  principle  of  the  heat 
pump  unit  is  similar  to  that  of  the  conventional  air  conditioning 
unit,  generally,  based  on  the  law  of  mass  and  energy  conservation, 
researchers  establish  mathematical  models  of  four  main  compo¬ 
nents  of  the  typical  heat  pump  unit,  including  the  compressor, 
condenser,  expansion  valve  and  evaporator  and  then  combine 
these  models  to  simulate  the  operation  performance  of  the  heat 
pump  unit  via  the  numerical  method.  Zhou  et  al.  [59]  proposed  a 
dynamic  lumped  parameter  model  to  describe  the  operating 
characteristics,  and  mass  and  energy  conservation  of  the  refrig¬ 
eration  equipment.  The  parameters  of  the  refrigeration  unit  in  any 
operating  mode  are  quickly  simulated  based  on  a  few  structural 
parameters  and  performance  parameters  of  the  working  condi¬ 
tion.  Yang  et  al.  [60]  developed  a  characteristic  parameter 
optimization  prediction  model,  which  is  used  to  analyze  the 
operation  performance  of  the  water-to-water  heat  pump  with 
eight  undetermined  characteristic  parameters.  Wang  [61] 
regarded  the  heating  output  of  the  heat  pump  unit  as  the  power 
function  of  the  flow  temperature  and  flow  rate  of  the  source  side 
and  load  side.  By  employing  the  neural  network  analyzing 
method,  the  operation  feature  curve  of  the  heat  pump  unit  is 
fitted  to  analyze  the  variation  of  its  heating  output  and  cooling 
output.  The  operation  feature  curve  also  illustrates  the  coupling 
relationship  between  the  ground  source  heat  pump  unit  and  the 
GHE  under  annual  dynamic  load.  Based  on  the  simulation  model 
of  the  GSHP  system  and  the  building  load  in  a  whole  year,  Yang 
et  al.  [62]  defined  the  annual  performance  factor  (APF)  of  the 
water-to-water  heat  pump  and  used  it  to  study  the  optimal 
design  and  matching  mode  between  different  components  for 
the  on/off  control  scheme  of  the  compressor. 

Some  researchers  adopted  a  more  convenient  calculation 
method  that  the  performance  of  the  heat  pump  unit  only  depends 
on  the  inlet  fluid  temperature.  Therefore,  the  performance  simu¬ 
lation  of  the  total  GSHP  system  is  achieved  by  combining  the 
outlet  fluid  temperature  of  the  GHE  and  the  operation  perfor¬ 
mance  of  the  heat  pump  unit  and  the  GHE.  Literatures  [63-66] 
fitted  the  functions  of  the  performance  of  the  heat  pump  unit 


and  the  circulating  fluid  outlet  temperature  of  the  GHE  according 
to  the  performance  data  of  the  heat  pump  unit  given  by 
manufacturers. 

At  present,  there  are  limited  studies  focusing  on  the  dynamic 
coupled  model  of  the  GHE  and  the  ground  heat  pump  unit  in 
China.  Wang  [67]  adopted  the  mechanism  modeling  method  to 
establish  a  dynamic  mathematical  model  combining  the  heat 
pump  unit  and  the  load  of  the  building.  Besides,  the  fuzzy  control 
method  is  applied  in  the  system.  The  dynamic  operating  simula¬ 
tion  is  carried  out  with  the  help  of  Matlab/Simulink.  Qu  [68]  built 
a  dynamic  simulation  model  by  incorporating  the  ground  heat 
pump  model,  GHE  model  and  water  loop  model  on  the  basis  of 
the  energy  and  mass  conservation.  The  simulation  model  can 
predict  the  operating  parameters  and  energy  consumption  profile 
of  the  GSHP  system  at  any  time. 

Zhao  et  al.  [69]  established  a  prediction  model  whose  key 
parameters  are  the  heat  exchange  rate  and  heat  transfer  coeffi¬ 
cient  by  using  the  grey  theory.  The  operational  conditions  of  the 
unbalance  of  the  load  between  winter  and  summer  are  divided 
into  three  working  phases:  transitional  phase,  steady  phase  and 
declining  phase.  Grey  prediction  can  provide  a  quantitative  basis 
for  the  stable  operation  and  the  operation  optimization  of  the 
GSHP  system.  At  the  mean  time,  experimental  investigation  [70] 
was  conducted  to  analyze  the  long-term  operation  performance 
of  the  GSHP.  The  experimental  results  showed  that  the  COP  of 
GSHP  system  increased  in  winter  and  decreased  in  summer.  Ma 
et  al.  [71]  detected  the  operation  of  an  actual  GSHP  project  for 
three  consecutive  quarters,  analyzing  the  reasons  and  influence 
factors  of  the  ground  heat  unbalance  and  proposing  advices  for 
the  optimization. 

It  is  known  that  the  soil  is  an  infinite  heat  carrier.  Centered  on 
the  GHE,  heat  is  released/extracted  to/from  surrounding  ground. 
The  load  accumulation  around  the  GHE  is  weakening  gradually 
and  the  capacity  of  heat  transfer  can  recover.  It  can  be  possible  if 
the  actual  project  is  intermittent  operation.  When  the  system  is 
closed  down,  the  capacity  of  heat  transfer  can  be  recovered  and 
the  performance  of  the  system  can  be  improved.  Many  research¬ 
ers  investigated  the  intermittent  operation  of  the  GSHP  system. 
Yang  et  al.  [72]  found  that  the  intermittent  operation  can  increase 
the  operating  efficiency  of  the  GSHP  system  by  the  experimental 
investigation.  By  controlling  the  operating  patterns  of  the  heat 
pump  unit,  Liu  et  al.  [73]  studied  the  changing  trend  of  the  water 
temperature  in  both  the  continuous  operating  mode  and  inter¬ 
mittent  operating  mode  in  order  to  find  out  the  best  approach  for 
optimizing  the  system  operation.  Wang  et  al.  [74]  indicated  that 
the  performance  of  the  GSHP  system  was  influenced  by  the  load 
of  the  building.  He  discussed  the  feasibility  of  the  intermittent 
operation  of  the  GSHP  system  by  considering  the  load  characteristics 
of  the  building.  Zhao  et  al.  [75]  carried  out  long-term  numerical 
simulation  on  the  heat  transfer  characteristics  of  the  soil  around 
the  foundation  pile  GHEs  covering  a  large  land  area,  presented  the 
concept  of  heat  barrier,  and  advised  that  the  system  analysis 
should  be  implemented  from  the  aspects  of  energy  equilibrium 
and  heat  barrier-to  ensure  the  whole  stability  and  reliability  of 
the  system.  Yu  et  al.  [76]  measured  the  actual  performance  of  the 
GHE  with  the  dissimilar  borehole  and  GHE  structures  and  found 
that  the  intermittent  operation  mode  can  maintain  a  better 
performance  of  GHEs  at  an  improved  value  of  33.9%  compared 
with  the  continuous  operation  condition.  Gao  et  al.  [77-79] 
analyzed  the  recovery  of  the  heat  transfer  capacity  of  the  soil 
when  the  GSHP  is  in  intermittent  process.  They  employed  the 
cylindrical  model  with  a  constant  heat  flux  to  calculate  the  soil 
temperature  recovery  during  the  intermittent  operation,  and 
performed  experiments  to  compare  with  the  ground  temperature 
changing  trend  in  the  intermittent  operation  and  continuous 
operation.  They  found  that  with  reasonable  manual  intermittent 
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control  based  on  the  actual  operational  conditions,  the  soil  heat 
transfer  capacity  can  be  increased.  Cui  et  al.  [80]  also  utilized  the 
superposition  principle  and  analytical  solution  to  simulate  the 
temperature  distribution  of  the  soil  around  the  GHE  in  the 
intermittent  operation.  Liu  et  al.  [81]  established  a  two-dimen¬ 
sional  heat  transfer  model  in  the  direction  of  radius  and  length, 
which  was  applied  to  calculate  the  recovery  condition  of  ground 
temperature  in  the  transitional  season  (spring).  Based  on  the 
layout  of  the  multiple  pipes  and  the  operating  conditions  of  an 
actual  GSHP  project,  Liu  et  al.  [82]  numerically  simulated  a  soil 
temperature  field  where  the  GHE  has  operated  for  5  years  in 
typical  area  and  under  typical  operational  conditions,  and  ana¬ 
lyzed  the  simulation  results  in  detail.  Yang  et  al.  [83]  presented 
the  soil  temperature  recovery  ratio  to  estimate  the  recovery 
status  of  the  soil  temperature. 

3.2.  Simulation  of  HGSHP 

To  deal  with  the  performance  decline  of  the  GSHP  system 
which  is  caused  by  the  unbalanced  building  load  between  winter 
and  summer,  the  HGSHP  with  assistant  cooling  devices  or  heating 
devices  has  appeared.  At  present,  most  studies  mainly  focus  on 
the  computational  analysis  of  the  operation  mode,  optimization 
and  performance  of  the  system.  Wang  et  al.  [84,85]  simulated  and 
evaluated  the  operation  performance  of  the  assistant  cooling- 
source  heat  pump  system  in  the  hot  summer  and  cold  winter 
region,  and  adopted  the  corrected  design  method  to  maximize  its 
performance  and  indicated  that  considering  the  initial  investment 
and  operating  cost,  the  assistant  cooling-source  heat  pump 
system  is  superior  to  the  GSHP  system  without  assistant  cooling. 
Yang  et  al.  [86,87]  developed  a  mathematic  model  for  the  solar- 
earth  source  heat  pump  system  (SESHPs)  on  the  basis  of  the 
component  models,  and  simulated  numerically  the  performance  of 
alternative  SESHPs.  Ma  et  al.  [88-93]  presented  the  integrated 
frozen  soil  cool  storage-GSHP  system.  Through  establishing  the 
mathematic  model  of  the  soil  charge  and  discharge,  the  operation 
performance  of  the  integrated  system  is  simulated  numerically  and 
the  energy  loss  through  the  soil  charge  and  discharge  is  analyzed.  In 
addition,  literatures  [94-96]  have  investigated  the  control  strate¬ 
gies,  design  procedures  and  optimization  designs  of  different  kinds 
of  HGSHP  systems  and  provided  general  methodologies  and  advices 
for  optimal  design  and  operation  of  such  HGSHP  systems. 

4.  Discussions  and  conclusions 

During  the  recent  few  decades,  great  deals  of  GSHP  systems 
have  been  widely  used  in  China.  In  this  paper,  most  heat  transfer 
simulation  models  currently  available  for  vertical  GHEs  have  been 
described  in  details;  the  experimental  investigation,  as  one  of  the 
important  methods  for  investigating  the  performance  of  the  GSHP 
system,  has  been  reviewed;  and  simulation  models  aiming  to 
evaluate  and  optimize  the  energy  performance  of  the  GSHP 
system  have  been  referred  as  well.  Finally,  this  paper  focuses  on 
the  simulation  of  the  HGSHP  system  with  assistant  cooling 
devices  or  heating  devices  since  the  HGSHP  system  is  the 
tendency  of  the  GSHP. 

It  can  be  seen  from  the  overview  that  Chinese  researchers  have 
obtained  the  following  achievements.  First,  the  traditional  infinite 
line  source  model  ignored  the  influence  of  the  finite  length  of 
borehole  and  the  boundary  condition  of  the  ground  surface.  On  the 
basis  of  the  infinite  line  source  model,  Chinese  researchers  utilized 
the  Green  function  method  to  obtain  the  temperature  response  of 
the  finite  line  heat  source  in  the  semi-infinite  medium.  At  the 
same  time,  the  equivalent  heat  transfer  resistances  of  the  single 
U-tube  GHE  and  double  U-tube  GHE  were  calculated,  and  the 


quasi  three-dimensional  heat  transfer  model  of  the  vertical  U-tube 
GHE  was  derived  by  using  the  superposition  principle.  And  the 
integral  average  temperature  along  the  length  direction  of  the 
borehole  wall  was  proposed  to  be  adopted  as  the  representative 
temperature  of  the  borehole  wall.  For  the  problem  of  the  ground- 
water  advection,  the  advection  of  the  porous  medium  was 
regarded  as  the  moving  heat  source  problem  and  the  analytical 
solution  of  the  temperature  response  caused  by  the  line  heat 
source  in  the  infinite  medium  with  advection  was  resolved. 

Second,  Chinese  researchers  proposed  computational  models 
by  combining  the  analytical  solution  and  numerical  solution  in 
order  to  give  considerations  to  both  the  computation  speed  and 
accuracy.  The  computational  domain  is  divided  into  two  parts 
based  on  the  boundary  of  the  borehole  wall.  The  computational 
method  is  the  steady  analytic  solution  for  the  borehole  inner  and 
the  finite  volume  method  for  the  soil.  The  simulation  models  of 
the  two  parts  are  coupled  based  on  the  borehole  wall  temperature 
and  the  heat  flux  through  the  borehole  wall. 

Although  a  great  number  of  researches  have  focused  on  the 
development  and  application  of  the  GSHP  system,  there  are  still 
a  few  aspects  that  need  further  investigation  to  broaden  and 
strengthen  the  applicability  of  the  GSHP  technology: 

(1)  The  vertical  GHEs  penetrate  through  different  structures  of 
geological  layers.  Both  the  multi-layered  soil  and  ground- 
water  have  great  influences  on  the  heat  transfer  performance 
of  vertical  GHEs.  However,  most  heat  transfer  models  of  GHEs 
cannot  take  into  the  effect  of  multi-layered  soil  and  ground- 
water.  Therefore,  the  studies  on  the  effects  of  these  two 
factors  should  be  strengthened  to  improve  current  heat 
transfer  models. 

(2)  In  actual  projects,  GHEs  are  always  made  up  of  multiple  pipes. 
It  is  difficult  to  accurately  simulate  the  instantaneous  perfor¬ 
mance  of  larger-scale  GHEs  because  the  simulating  calcula¬ 
tion  is  very  complex  and  time-consuming.  There  are  few 
reports  aiming  to  optimize  the  operation  strategies  of  multi¬ 
ple  pipes  in  partial  load.  It  is  worth  to  carry  out  investigations 
on  this  issue  in  the  future. 

(3)  During  the  operation  of  the  GSHP  system,  the  building  load, 
the  heat  pump  unit  and  the  GHE  influence  each  other.  So,  it  is 
necessary  to  set  up  the  coupled  heat  transfer  model  of  the 
ground  source  heat  pump  system  to  obtain  more  accurate 
dynamic  performance  of  the  system. 
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